BCG analyzed

install packages

#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")

Load libraries

In this step, we are loading the libraries which are required for analysis and plotting of the data.

library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(limma)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)

set the colors of the cell types

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
    'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")

load data for mouse 1 and mouse 2

# Mouse 1
bcg_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/bcg_sample_round_one/outs/"

list.files(bcg_mouse_1)

##  [1] "analysis"                      "cloupe.cloupe"                
##  [3] "filtered_feature_bc_matrix"    "filtered_feature_bc_matrix.h5"
##  [5] "metrics_summary.csv"           "molecule_info.h5"             
##  [7] "possorted_genome_bam.bam"      "possorted_genome_bam.bam.bai" 
##  [9] "probe_set.csv"                 "raw_feature_bc_matrix"        
## [11] "raw_feature_bc_matrix.h5"      "spatial"                      
## [13] "spatial_enrichment.csv"        "web_summary.html"

bcg_mouse_1_data <- Load10X_Spatial(
  bcg_mouse_1,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "bcg",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

# Mouse 2
bcg_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/bcg_sample_round_two/outs/"

list.files(bcg_mouse_2)

##  [1] "analysis"                      "bcg_only_web_summary.html"    
##  [3] "cloupe.cloupe"                 "filtered_feature_bc_matrix"   
##  [5] "filtered_feature_bc_matrix.h5" "metrics_summary.csv"          
##  [7] "molecule_info.h5"              "possorted_genome_bam.bam"     
##  [9] "possorted_genome_bam.bam.bai"  "probe_set.csv"                
## [11] "raw_feature_bc_matrix"         "raw_feature_bc_matrix.h5"     
## [13] "spatial"                       "spatial_enrichment.csv"

bcg_mouse_2_data <- Load10X_Spatial(
 bcg_mouse_2,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "bcg",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

data preprocessing

# Mouse 1
bcg_mouse_1_data <- SCTransform(bcg_mouse_1_data, assay = "spatial", verbose = FALSE)

# Mouse 2
bcg_mouse_2_data <- SCTransform(bcg_mouse_2_data, assay = "spatial", verbose = FALSE)

Visualization of lymphoid follicle markers

Mouse 1

p1 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(bcg_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

Mouse 2

p1 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(bcg_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

dimensionality reduction and clustering

# Mouse 1
bcg_mouse_1_analyzed <- RunPCA(bcg_mouse_1_data, assay = "SCT", verbose = FALSE)
bcg_mouse_1_analyzed <- FindNeighbors(bcg_mouse_1_analyzed, reduction = "pca", dims = 1:5)
bcg_mouse_1_analyzed <- FindClusters(bcg_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
bcg_mouse_1_analyzed <- RunUMAP(bcg_mouse_1_analyzed, reduction = "pca", dims = 1:5)

# Mouse 2
bcg_mouse_2_analyzed <- RunPCA(bcg_mouse_2_data, assay = "SCT", verbose = FALSE)
bcg_mouse_2_analyzed <- FindNeighbors(bcg_mouse_2_analyzed, reduction = "pca", dims = 1:5)
bcg_mouse_2_analyzed <- FindClusters(bcg_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
bcg_mouse_2_analyzed <- RunUMAP(bcg_mouse_2_analyzed, reduction = "pca", dims = 1:5)

labelling clusters using SingleR pipeline

# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)

# Mouse 1
bcg_mouse_1_analyzed.sce <- as.SingleCellExperiment(bcg_mouse_1_analyzed)

# Mouse 2
bcg_mouse_2_analyzed.sce <- as.SingleCellExperiment(bcg_mouse_2_analyzed)

#create reference data

ref.data <- celldex::ImmGenData()

# run singleR pipeline to find the cell types

## Mouse 1
cell_types_mouse_1 <- SingleR(test=bcg_mouse_1_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

## Mouse 2
cell_types_mouse_2 <- SingleR(test=bcg_mouse_2_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

Check if the cell types are labelled with high accuracy

# view cell types

# Mouse 1
table(cell_types_mouse_1$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##                48                43              1244               264 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               437              1052                20                 2 
##     Stromal cells 
##               208

#checking cell types Mouse 1

plotScoreHeatmap(cell_types_mouse_1)

plotDeltaDistribution(cell_types_mouse_1)

summary(is.na(cell_types_mouse_1$pruned.labels))

##    Mode   FALSE 
## logical    3318

# Mouse 2
table(cell_types_mouse_2$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##               130                62              1027               175 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               335               707                 5                 1 
##     Stromal cells 
##               258

#checking cell types Mouse 2

plotScoreHeatmap(cell_types_mouse_2)

plotDeltaDistribution(cell_types_mouse_2)

summary(is.na(cell_types_mouse_2$pruned.labels))

##    Mode   FALSE    TRUE 
## logical    2694       6

merging metadata of cell types to seurat object

# Mouse 1
bcg_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels

# Mouse 2
bcg_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels

Plot umap and spatial plot to localize each cell type

Mouse 1

p1 <- DimPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Mouse 2

p1 <- DimPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Plot violin plot with top markers associated with B cells and lymphoid follicles

Mouse 1

VlnPlot(bcg_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

Mouse 2

VlnPlot(bcg_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

find markers in each cluster

Mouse 1

B_cell_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)

##                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cr2       8.320332e-60  0.9384812 0.604 0.052 1.291315e-55
## Ms4a1     1.615581e-56  0.8401004 0.604 0.055 2.507382e-52
## Ighd      4.694334e-53  1.3861847 0.812 0.120 7.285606e-49
## Bank1     8.505742e-52  0.7844861 0.646 0.068 1.320091e-47
## Cd19      4.140845e-51  1.4892818 0.917 0.167 6.426591e-47
## Il9r      3.138481e-50  0.4293062 0.354 0.020 4.870923e-46
## Blnk      3.244900e-45  1.0068607 0.812 0.132 5.036085e-41
## Tnfrsf13c 1.132859e-43  1.1269492 0.750 0.119 1.758198e-39
## Spib      2.486053e-39  1.2094113 0.708 0.117 3.858354e-35
## Pla2g2d   2.630705e-38  0.7144449 0.479 0.050 4.082853e-34

write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/B_cell_mouse_1.csv")

DC_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Fam241a 1.741273e-09  0.2699921 0.279 0.058 2.702456e-05
## Slfn9   6.331009e-09  0.2754009 0.256 0.053 9.825725e-05
## Fmo2    9.853839e-09 -0.9743768 0.349 0.690 1.529316e-04
## Epas1   1.523917e-08 -0.8608185 0.605 0.896 2.365120e-04
## Plxnc1  3.464226e-08  0.4707331 0.535 0.199 5.376479e-04
## Igfbp2  7.898749e-08 -0.8938734 0.558 0.867 1.225886e-03
## Slpi    4.991795e-07  0.4812667 0.953 0.868 7.747265e-03
## Aqp1    9.984770e-07 -0.7826018 0.233 0.595 1.549636e-02
## Ltb     1.031783e-06  0.3426673 0.465 0.175 1.601327e-02
## Sftpc   1.082770e-06 -0.2978901 1.000 0.999 1.680459e-02

write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/DC_mouse_1.csv")

Stromal_cells_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igkc    2.549004e-17 -0.5568178 1.000 0.999 3.956055e-13
## Igha    1.333650e-16 -0.5432519 1.000 0.997 2.069825e-12
## Gsn     1.929346e-15  0.4734754 0.933 0.764 2.994344e-11
## Jchain  5.250575e-15 -0.6192559 0.654 0.849 8.148893e-11
## Tnfaip2 2.023024e-13 -0.4416083 0.933 0.964 3.139733e-09
## Epas1   3.849126e-12  0.3735630 0.995 0.885 5.973844e-08
## Inmt    4.073284e-12  0.4501445 0.981 0.841 6.321736e-08
## Ighj4   1.951249e-11 -0.5177593 0.630 0.801 3.028338e-07
## Timp3   2.877549e-11  0.3253112 0.995 0.967 4.465955e-07
## Cldn5   2.079806e-10  0.3490426 0.952 0.835 3.227860e-06

write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Stromal_cells_mouse_1.csv")


Endothelial_cells_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igfbp2  2.960434e-64  0.6943825 0.940 0.817 4.594593e-60
## Epas1   4.555557e-63  0.5417082 0.957 0.853 7.070224e-59
## Cldn5   7.254164e-62  0.5775664 0.934 0.788 1.125846e-57
## Tns1    3.981420e-47  0.4378691 0.965 0.883 6.179164e-43
## Ace     1.107643e-46  0.4589540 0.747 0.528 1.719061e-42
## Spock2  2.544447e-44  0.4939606 0.722 0.540 3.948982e-40
## Pakap.1 6.648450e-44  0.4067602 0.920 0.843 1.031840e-39
## Rasip1  6.994426e-43  0.4084629 0.716 0.498 1.085535e-38
## Timp3   9.312203e-43  0.3806280 0.994 0.953 1.445254e-38
## Clec14a 4.269361e-42  0.3756253 0.589 0.358 6.626049e-38

write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Endothelial_cells_mouse_1.csv")

Epithelial_cell_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Inmt    3.735506e-25 -0.8116861 0.678 0.865 5.797505e-21
## Tns1    1.355276e-21 -0.5743788 0.818 0.922 2.103388e-17
## Igha    2.119608e-20  0.3895227 1.000 0.997 3.289632e-16
## Igfbp2  2.733816e-19 -0.7586417 0.712 0.876 4.242883e-15
## Cldn5   2.164530e-18 -0.6047003 0.708 0.854 3.359351e-14
## Epas1   4.393914e-18 -0.5764613 0.777 0.902 6.819354e-14
## Timp3   1.315169e-17 -0.4704262 0.928 0.972 2.041142e-13
## Pakap.1 3.745781e-17 -0.4711075 0.777 0.880 5.813452e-13
## Igkc    4.927450e-17  0.3794999 1.000 0.999 7.647402e-13
## Cd93    7.279343e-17 -0.4278941 0.178 0.446 1.129754e-12

write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Epithelial_cell_mouse_1.csv")

Fibroblasts_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Inmt    5.509154e-26  0.4909881 0.966 0.832 8.550208e-22
## Gsn     1.321252e-25  0.4882381 0.906 0.754 2.050583e-21
## Ltbp4   1.503307e-17  0.3913844 0.654 0.474 2.333132e-13
## Clec3b  1.204074e-16  0.4400577 0.613 0.437 1.868723e-12
## Nbl1    2.244220e-15  0.3431696 0.883 0.774 3.483029e-11
## Fmo2    1.080310e-14  0.3284504 0.808 0.667 1.676641e-10
## Rarres2 2.062960e-14  0.4293750 0.751 0.625 3.201714e-10
## Timp3   1.157632e-13  0.2623184 0.993 0.965 1.796645e-09
## Gbp2b   2.138363e-13 -0.2928538 0.924 0.953 3.318739e-09
## Ogn     5.905678e-13  0.2581432 0.384 0.234 9.165613e-09

  write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Fibroblasts_mouse_1.csv")
  
Macrophages_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igfbp2  2.890373e-66 -0.7985205 0.787 0.898 4.485859e-62
## Epas1   1.058095e-59 -0.5990300 0.828 0.922 1.642164e-55
## Igha    7.269746e-58  0.4226417 0.997 0.997 1.128265e-53
## Inmt    1.475707e-55 -0.7200217 0.775 0.885 2.290297e-51
## Cldn5   5.228715e-52 -0.6074341 0.761 0.880 8.114965e-48
## Igkc    3.122946e-49  0.3467899 0.998 0.999 4.846812e-45
## Ctss    3.077046e-48  0.4188831 0.998 0.981 4.775575e-44
## Timp3   8.672736e-48 -0.4636797 0.952 0.976 1.346009e-43
## Tnfaip2 1.139208e-45  0.3668298 0.993 0.948 1.768051e-41
## Tns1    3.488936e-43 -0.4879362 0.872 0.934 5.414828e-39

write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophages_mouse_1.csv")

Monocytes_mouse_1 <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Monocytes_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igfbp2  2.890373e-66 -0.7985205 0.787 0.898 4.485859e-62
## Epas1   1.058095e-59 -0.5990300 0.828 0.922 1.642164e-55
## Igha    7.269746e-58  0.4226417 0.997 0.997 1.128265e-53
## Inmt    1.475707e-55 -0.7200217 0.775 0.885 2.290297e-51
## Cldn5   5.228715e-52 -0.6074341 0.761 0.880 8.114965e-48
## Igkc    3.122946e-49  0.3467899 0.998 0.999 4.846812e-45
## Ctss    3.077046e-48  0.4188831 0.998 0.981 4.775575e-44
## Timp3   8.672736e-48 -0.4636797 0.952 0.976 1.346009e-43
## Tnfaip2 1.139208e-45  0.3668298 0.993 0.948 1.768051e-41
## Tns1    3.488936e-43 -0.4879362 0.872 0.934 5.414828e-39

write.csv(Monocytes_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Monocytes_mouse_1.csv")

Mouse 2

B_cell_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)

##                   p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cxcr5     1.795952e-103  0.8833483 0.638 0.072 2.455426e-99
## Tnfrsf13c  1.976759e-98  1.6255464 0.877 0.182 2.702624e-94
## Cr2        2.568990e-98  1.1656010 0.623 0.075 3.512323e-94
## Cd19       2.894342e-97  1.6119812 0.885 0.195 3.957144e-93
## Iglc3      2.116963e-92  1.2005717 0.862 0.178 2.894312e-88
## Ighd       3.359851e-88  1.4803253 0.808 0.165 4.593588e-84
## Fcrla      1.286344e-85  0.9721240 0.677 0.105 1.758689e-81
## Bank1      1.955197e-84  1.0070053 0.692 0.112 2.673146e-80
## Pax5       2.727178e-83  0.9189419 0.662 0.102 3.728597e-79
## Ly6d       1.258899e-81  0.5611298 0.515 0.056 1.721167e-77

write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/B_cell_mouse_2.csv")

DC_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Itk    1.126374e-24  0.6518767 0.806 0.240 1.539978e-20
## Sftpc  2.294687e-22 -1.2702568 1.000 1.000 3.137296e-18
## Pak1   2.942115e-20  0.5473896 0.726 0.218 4.022460e-16
## Clic5  4.911070e-20 -1.2458544 0.694 0.956 6.714415e-16
## C3     5.223048e-20  0.8231177 1.000 0.997 7.140951e-16
## Cd3d   5.577503e-20  0.8201085 0.839 0.343 7.625563e-16
## Ccr7   1.958530e-19  0.6471925 0.629 0.187 2.677702e-15
## Ager   5.354944e-19 -1.1661539 0.903 0.973 7.321280e-15
## H2-Eb1 7.763557e-19  0.6541973 1.000 0.999 1.061433e-14
## Vegfa  1.098108e-18 -1.1992835 0.629 0.932 1.501333e-14

write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/DC_mouse_2.csv")

Stromal_cells_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Gsn     1.386259e-28  0.6018293 0.996 0.905 1.895294e-24
## Tnfaip2 4.346716e-23 -0.7090419 0.953 0.984 5.942831e-19
## Psap    1.665689e-21 -0.4323633 1.000 1.000 2.277331e-17
## Myl9    1.259329e-19  0.4828020 0.477 0.239 1.721755e-15
## H2-DMb2 1.679277e-19 -0.6117161 0.946 0.977 2.295907e-15
## Ctss    4.810568e-19 -0.5572010 0.992 0.995 6.577009e-15
## Ly6e    5.775611e-19 -0.4271559 0.984 0.998 7.896415e-15
## B2m     1.406355e-18 -0.3576475 0.996 0.998 1.922769e-14
## Cyp2f2  7.401773e-18  0.8529229 1.000 0.941 1.011970e-13
## Myh11   3.177795e-17  0.5537961 0.473 0.249 4.344681e-13

write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Stromal_cells_mouse_2.csv")


Endothelial_cells_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)

##                 p_val avg_log2FC pct.1 pct.2     p_val_adj
## Cldn5   8.891607e-143  0.9344905 0.985 0.738 1.215661e-138
## Igfbp2  4.848456e-140  1.0061615 0.983 0.708 6.628810e-136
## Tns1    1.456044e-130  0.7875292 1.000 0.882 1.990704e-126
## Epas1   2.648369e-124  0.7706319 0.997 0.854 3.620850e-120
## Inmt    2.626647e-108  0.8388433 0.979 0.682 3.591152e-104
## Timp3   2.201202e-102  0.6322175 0.998 0.929  3.009484e-98
## Clic5   9.961688e-101  0.5969119 0.999 0.921  1.361962e-96
## Pakap.1  1.276464e-99  0.6482757 0.990 0.892  1.745182e-95
## Spock2   1.998597e-97  0.7174474 0.937 0.628  2.732482e-93
## Acer2    2.689486e-93  0.6562824 0.951 0.763  3.677065e-89

write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Endothelial_cells_mouse_2.csv")

Epithelial_cell_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)

##                      p_val avg_log2FC pct.1 pct.2    p_val_adj
## Igfbp2        4.276103e-42 -1.7397435 0.474 0.836 5.846287e-38
## Inmt          4.796397e-41 -1.6731033 0.474 0.817 6.557634e-37
## Spint2        3.086740e-40  0.7716912 1.000 0.972 4.220190e-36
## Ccdc96        9.863237e-40  0.4891312 0.326 0.059 1.348502e-35
## Wfdc2         3.877575e-38  1.3042103 0.943 0.806 5.301421e-34
## Tns1          2.295413e-37 -1.1410230 0.783 0.937 3.138288e-33
## Cdhr4         3.244769e-37  0.7131494 0.400 0.096 4.436249e-33
## Erich3        6.609781e-37  0.4512833 0.320 0.060 9.036892e-33
## 1700001C02Rik 1.067681e-36  0.3713457 0.303 0.053 1.459733e-32
## Cldn5         1.185186e-36 -1.3527028 0.583 0.849 1.620386e-32

write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Epithelial_cell_mouse_2.csv")

Fibroblasts_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cst3    4.019251e-33  0.4103062 1.000 1.000 5.495120e-29
## Gsn     6.675412e-26  0.5246076 0.976 0.905 9.126623e-22
## Clec3b  5.877289e-24  0.5597167 0.887 0.725 8.035430e-20
## Mgp     1.809036e-23  0.3799454 1.000 0.992 2.473314e-19
## Ltbp4   9.015487e-22  0.5156064 0.866 0.682 1.232597e-17
## Inmt    1.181286e-21  0.5020674 0.961 0.771 1.615054e-17
## Eln     6.146593e-17  0.4589563 0.815 0.668 8.403622e-13
## H2-DMb2 2.324752e-16 -0.4986663 0.967 0.975 3.178401e-12
## H2-Ab1  1.293291e-15 -0.3523212 0.988 0.991 1.768188e-11
## Ctss    1.964660e-15 -0.4255406 0.997 0.995 2.686083e-11

  write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Fibroblasts_mouse_2.csv")
  
Macrophages_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Ctss    1.040821e-95  0.7567875 1.000 0.993 1.423011e-91
## Psap    3.925294e-89  0.5558869 1.000 1.000 5.366662e-85
## Epas1   2.949012e-88 -0.8781662 0.827 0.937 4.031889e-84
## Inmt    6.472508e-85 -1.1013762 0.611 0.860 8.849213e-81
## Gbp2b   9.581003e-81  0.6528410 0.993 0.962 1.309915e-76
## Tnfaip2 2.300836e-79  0.7575139 0.997 0.976 3.145703e-75
## Cldn5   4.546197e-76 -0.9452875 0.699 0.879 6.215561e-72
## Ctsc    1.261996e-71  0.6238895 0.997 0.957 1.725401e-67
## Nos2    1.006443e-70  0.9034125 0.755 0.441 1.376009e-66
## Igfbp2  4.797133e-70 -1.0009458 0.666 0.865 6.558640e-66

write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophages_mouse_2.csv")

Monocytes_mouse_2 <- FindMarkers(bcg_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Monocytes", min.pct = 0.25)
  head(Monocytes_mouse_2, n = 10)

##                      p_val avg_log2FC pct.1 pct.2    p_val_adj
## Epha8         1.726519e-32  0.4779517   0.4 0.005 2.360496e-28
## Zcchc18       1.701735e-23  0.4752914   0.4 0.007 2.326611e-19
## Nxpe4         3.777130e-19  0.4726359   0.4 0.009 5.164092e-15
## Saa3          1.956962e-12  0.4641713   0.4 0.014 2.675559e-08
## Cacng8        3.698153e-11  0.4610098   0.4 0.016 5.056115e-07
## Tmem69        5.716732e-10  0.4589060   0.4 0.019 7.815916e-06
## Cd5           8.543622e-10  0.4583805   0.4 0.019 1.168084e-05
## Paip1         3.760272e-09  0.4557560   0.4 0.020 5.141044e-05
## Olfm2         4.327556e-09  0.4494765   0.4 0.020 5.916634e-05
## 1810058I24Rik 1.008572e-08  0.4536598   0.4 0.022 1.378920e-04

write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Monocytes_mouse_2.csv")

Macrophage_cell_deseq <- FindMarkers(bcg_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(Macrophage_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophage_cell_deseq.csv")

Cell-cell interaction using cell chat

bcg_cellchat <- createCellChat(object = bcg_mouse_2_analyzed, group.by = "ref.data")

## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information 
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells

CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
bcg_cellchat@DB <- CellChatDB.use


bcg_cellchat <- subsetData(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedGenes(bcg_cellchat)
bcg_cellchat <- identifyOverExpressedInteractions(bcg_cellchat)


bcg_cellchat <- computeCommunProb(bcg_cellchat)
bcg_cellchat <- filterCommunication(bcg_cellchat, min.cells = 10)

## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  Monocytes Neutrophils

bcg_cellchat <- computeCommunProbPathway(bcg_cellchat)
bcg_cellchat <- aggregateNet(bcg_cellchat)

groupSize <- as.numeric(table(bcg_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(bcg_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(bcg_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

# See interaction of each cell type with the another

mat1 <- bcg_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}

##Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
bcg_cellchat@netP$pathways

##  [1] "COLLAGEN"    "APP"         "GALECTIN"    "MHC-II"      "FN1"        
##  [6] "THBS"        "MIF"         "LAMININ"     "GAS"         "VEGF"       
## [11] "COMPLEMENT"  "SEMA4"       "SEMA3"       "JAM"         "ICAM"       
## [16] "CCL"         "ITGAL-ITGB2" "SPP1"        "CD22"        "CD45"       
## [21] "PDGF"        "TGFb"        "CXCL"        "GRN"         "CHEMERIN"   
## [26] "PARs"        "TENASCIN"    "PECAM1"      "WNT"         "TNF"        
## [31] "NOTCH"       "CD52"        "PD-L1"       "ESAM"        "SEMA7"      
## [36] "CDH"         "EPHA"        "IL16"        "THY1"        "ncWNT"      
## [41] "IGF"         "APRIL"       "ICOS"        "TWEAK"       "PERIOSTIN"  
## [46] "CDH5"        "PDL2"

# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
?netVisual_aggregate(bcg_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(bcg_cellchat, signaling = pathways.show, layout = "chord")

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(bcg_cellchat, signaling = pathways.show, color.heatmap = "Reds")

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(bcg_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(bcg_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(bcg_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(bcg_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(bcg_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)

Compute the network centrality scores

bcg_cellchat <- netAnalysis_computeCentrality(bcg_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(bcg_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(bcg_cellchat, pattern = "incoming")
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(bcg_cellchat, pattern = "outgoing")

nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(bcg_cellchat, pattern = "outgoing")

# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "outgoing")

Identify and visualize communication patterns (incoming signal)

selectK(bcg_cellchat, pattern = "incoming")

nPatterns = 3
bcg_cellchat <- identifyCommunicationPatterns(bcg_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(bcg_cellchat, pattern = "incoming")

# dot plot
netAnalysis_dot(bcg_cellchat, pattern = "incoming")

Gene set enrichment

organism = org.Mm.eg.db

prepare input data for gene set enrichment

# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/bcg/Macrophage_cell_deseq.csv", header=TRUE)

# we want the log2 fold change 
original_gene_list <- df$avg_log2FC

# name the vector
names(original_gene_list) <- df$X

# omit any NA values 
gene_list<-na.omit(original_gene_list)

# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)

Gene set enrichment

#gene set enrichment function for B cell cluster

gse <- gseGO(geneList=gene_list, 
             ont ="ALL", 
             keyType = "SYMBOL", 
             minGSSize = 3, 
             maxGSSize = 800, 
             pvalueCutoff = 0.05, 
             verbose = TRUE, 
             OrgDb = organism, 
             pAdjustMethod = "none",
             eps = 0.0)

dotplot(gse, showCategory=10, split=".sign") + facet_grid(.~.sign)